Accurate simulations of the dynamical bar-mode instability in full General Relativity 
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We present accurate simulations of the dynamical bar-mode instability in full General Relativity focussing 
on two aspects which have not been investigated in detail in the past. Namely, on the persistence of the bar 
deformation once the instability has reached its saturation and on the precise determination of the threshold for 
the onset of the instability in terms of the parameter l3 — T/\W\. We find that generic nonlinear mode-coupling 
effects appear during the development of the instability and these can severely limit the persistence of the bar 
deformation and eventually suppress the instability. In addition, we observe the dynamics of the instability 
to be strongly influenced by the value f3 and on its separation from the critical value /3c marking the onset 
of the instability. We discuss the impact these results have on the detection of gravitational waves from this 
process and provide evidence that the classical perturbative analysis of the bar-mode instability for Newtonian 
and incompressible Maclaurin spheroids remains qualitatively valid and accurate also in full General Relativity. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.40.Dg, 95.30.Lz, 97.60.Jd 95.30.Sf 



I. INTRODUCTION 

It is well known that rotating neutron stars are subject 
to non-axisymmetric instabilities for non-radial axial modes 
with azimuthal dependence e"""^ (with m = 1,2,...) when 
the instability parameter /3 = T'/lW^I ('-e. the ratio between 
the rotational kinetic energy T and the gravitational binding 
energy W) exceeds a critical value /3c- 

An exact and perturbative treatment of these instabilities 
exists only for incompressible self-gravitating fluids in New- 
tonian gravity (see refs. 01 01) and this predicts that a dynam- 
ical instability should arise for the "bar-mode" (i.e. the one 
with m=2) when (3 > /3c ^ 0.2738. On the other hand, the 
accurate study of the dynamical bar-mode instability, of its 
nonlinear evolution and of the determination of the threshold 
for the instability, demand the use of numerical simulations 
with the solution in three spatial dimensions (3D) of the fully 
nonlinear hydrodynamical equations coupled to the Einstein 
field equations. 

Despite these requirements, much of the literature on this 
process has so far been limited to a Newtonian or post- 
Newtonian (PN) description. While this represents an ap- 
proximation, these studies have provided important informa- 
tion on several aspects of the instability that could not have 
been investigated with perturbative techniques. In particu- 
lar, these numerical studies have shown that (3c depends very 
weakly on the stiffness of the equation of state (EOS) and that, 
once a bar has developed, the formation of spiral arms is im- 
portant for the redistribution of the angular momentum (see 
refs. 11111111,111111 [Mill). More recently, instead, it was 
shown that the threshold for the onset of the dynamical insta- 
bility can be smaller for stars with a high degree of differential 
rotation and a weak dependence on the EOS was confirmed in 
refs. 01 [3 [M [S [Hi- Finally, these Newtonian analy- 
ses have also provided the first evidence that an m=l-mode 
dynamical instability may play an important role for smaller 
values of the critical parameter This is also referred to as the 



"low-/3" instability OlTlllSIl and it will not be considered here. 

Only very recently it has become possible to perform simu- 
lations of the dynamical bar instability for old neutron stars 
in full General Relativity [19]. These studies have shown 
that within a fully general-relativistic framework the critical 
value for the onset of the instability is smaller than the New- 
tonian one (i.e. f3c — 0.24 — 0.25) and this behaviour was 
confirmed by PN calculations I20l l2lll which also suggested 
that (3c varies with the compactness M/R of the star. 

The bar-mode instability may take place in young neutron 
stars, either as the result of the accretion-induced collapse of 
a white dwarf ||22ll or in the collapse of a massive stellar core. 
Indeed, recent simulations investigating axisymmetric stellar- 
core collapse in full General Relativity [23] have pointed out 
that for sufficiently differentially rotating progenitors, it is, at 
least in principle, possible to obtain toroidal protoneutron-star 
cores with masses between 1.2 and 2.1 Mq which are unsta- 
ble against bar-mode deformations (it should be mentioned 
that it is still unclear how likely these high rotation-rates and 
strongly differential-rotation profiles actually are in nature). 
Besides pointing out this alternative interesting scenario, the 
work in ref. [23] has also suggested that more realistic EOS 
and neutrino cooling could enhance the process; this scenario 
has also been considered in the PN approximation in ref. ||24ll. 
An important common feature in all these investigations is that 
the development of the bar-mode instability was obtained in- 
troducing very strong ad hoc m=2-mode perturbations. 

The recent general-relativistic studies, together with their 
Newtonian and PN counterparts, have been very helpful in 
highlighting the main features of the instability. However, 
several fundamental questions remain unanswered. Most no- 
tably: i) What is the role of the initial perturbation? ii) What 
is the effect of the symmetry conditions often used in numeri- 
cal calculations? Hi) How are the dynamics influenced by the 
value of the parameter (3, especially when this is largely over- 
critical? Finally and most importantly: iv) How long does a 
bar survive, once fully developed? 
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Clearly, the last question has important implications for the 
possible observational relevance of the gravitational-wave sig- 
nal emitted through the bar-mode instability as the signal-to- 
noise ratio (SNR) can increase considerably in the case of a 
long-lived bar since the SNR grows as the square root of the 
number of the effective cycles of the signal available for de- 
tection. Earlier work on this subject basically suggested that 
once a bar was formed it would tend to be persistent on the 
radiation-reaction timescale. Indications and evidences in this 
direction were presented in a Newtonian framework in ref. f3\ 
as well as in a PN one in ref. |21]. In contrast, fully general- 
relativistic results, either using suitable symmetry boundary 
conditions [i.e. the so called 7r-symmetry boundary condi- 
tions) ifigll or not ll23ll . show a non-persistent bar 

In this work we try to find answers to these important open 
questions by exploring, in a systematic way, the bar-mode in- 
stability for a large number of initial stellar models. In doing 
this, we intend to go beyond the standard phenomenological 
discussion of the nonlinear dynamics of the instability often 
encountered in the literature. 

The main results of our analysis can be summarised as fol- 
lows: i) The initial perturbation (either in the form of an m=l 
mode or of an m=2 mode) can play a role in determining the 
duration of the bar-mode deformation, but not in determin- 
ing the growth time of the instability; the only exception to 
this is represented by models near the threshold; //) For mod- 
erately overcritical models (i.e. with /3 ~ f3c), the use of a 
TT-symmetry can radically change the dynamics and extend 
considerably the persistence of the bar; this ceases to be true 
for largely overcritical models (i.e. with f3 ^ /3c), for which 
even the artificial symmetries are not sufficient to provide a 
long-lived bar; The persistence of the bar is strongly de- 
pendent on the degree of overcriticality and is generically of 
the order of the dynamical timescale; iv) Generic nonlinear 
mode-coupling effects (especially between the m=l and the 
m=2 modes) appear during the development of the instability 
and these can severely limit the persistence of the bar defor- 
mation and eventually suppress the bar deformation; v) The 
dynamics of largely overcritical models are fully determined 
by the excess of rotational energy and the bar deformation 
is very rapidly suppressed through the conversion of kinetic 
energy into internal one. In addition, we have also assessed 
the accuracy of the classical Newtonian stability analysis of 
Maclaurin spheroids for incompressible self-gravitating flu- 
ids Overall, and despite having applied it to differentially 
rotating and relativistic models, we have found it to be sur- 
prisingly accurate in determining both the threshold for the 
instability and the complex eigenfrequencies for the unstable 
models. 

The paper is organized as follows. In Sec. II we give de- 
tails on the evolution methods used, while in Sec. Ill we dis- 
cuss the initial models and their properties. In Sec. IV we 
introduce the methodology used to analyse the numerical re- 
sults of the simulations, which are then discussed in Sec. V 
in terms of the general dynamics of the instability and of the 
general properties. In Sec. VI, we present the features of the 
instability that are specific to different treatments of the ini- 
tial conditions, while in Sec. VII we illustrate two different 



methods for the determination of Pc- Finally, in Sec. VIII we 
discuss the impact of our results on the emission of gravita- 
tional waves from the unstable models and present in Sec. IX 
our conclusions and the prospects of future research. 

We have used a space like signature (—,+,+,+), with 
Greek indices running from to 3, Latin indices from 1 to 
3 and the standard convention for the summation over re- 
peated indices. Furthermore, we indicate as (x,y,z) the Carte- 
sian coordinates and we define r ~ \/ + y'^ + z^, w = 
+ 6 ~ arctan(tu/z), cf) — arctan(y/x) for the ax- 
ial and spherical coordinates. Unless explicitly stated, all the 
quantities are expressed in the system of adimensional units 
in which c = G = Mq — 1. 



II. EVOLUTION OF FIELDS AND MATTER 

The code and the evolution method are the same as the ones 
used in Baiotti et al. |25', '26'] and therein described. For con- 
venience we report here the main properties and characteris- 
tics of the employed simulation method. We have used the 
general-relativistic hydrodynamics code Whisky, in which 
the hydrodynamics equations are written as finite differences 
on a Cartesian grid and solved using high-resolution shock- 
capturin g (H RSC) schemes (a first description of the code was 
given in ll26ll ). 



A. Evolution of Einstein equations 

The original ADM formulation casts the Einstein equations 
into a first-order (in time) quasi-linear iIztIi system of equa- 
tions. The dependent variables are the three-metric 7^ and 
the extrinsic curvature Kij, with first-order evolution equa- 
tions given by 



(2.1) 



(2.2) 

Here, a is the lapse function, (3i is the shift vector, denotes 
the CO variant derivative with respect to the three-metric 7^, 
Rij is the Ricci curvature of the three-metric, K = Kij 
is the trace of the extrinsic curvature, Sij is the projection 
of the stress-energy tensor onto the space-like hypersurfaces 
and S = 1^-' Sij (for a more detailed discussion, see ^28*]). 
In addition to the evolution equations, the Einstein equations 
also provide four constraint equations to be satisfied on each 
space-like hypersurface. These are the Hamiltonian constraint 
equation 



167rp^ 



, (2.3) 
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and the momentum constraint equations 

VjK'^ --/'^VjK -8nf ^0 . (2.4) 

In equations ( l2.1l i-( IZ4l i. p^dm ^^'^ J* energy density 

and the momentum density as measured by an observer mov- 
ing orthogonally to the space-like hypersurfaces. 

In particular, we use a conformal traceless reformulation of 
the above system of evolution equations, as first suggested by 
Nakamura, Oohara and Kojima | |29t] (NOK formulation), in 
which the evolved variables are the conformal factor (0), the 
trace of the extrinsic curvature (K), the conformal 3-metric 
{jij), the conformal traceless extrinsic curvature (A^j) and 
the conformal connection functions (F*), defined as 

= ilog(^) , (2.5) 

K = f^K,, , (2.6) 

% = e-"*7*j- , (2.7) 

= e-^^{K,j--i,,K) , (2.8) 

= 7;"^' . (2.9) 

The code used for evolving these quantities is the one de- 
veloped within the Cactus computational toolkit [30] and is 
designed to handle arbitrary shift and lapse conditions. In par- 
ticular, we have used hyperbolic if -driver slicing conditions 
of the form 

dta = -f{a)a^{K-Ko), (2.10) 

with f{a) > and Kq = K{t = 0). This is a generalization 
of many well known slicing conditions. For example, setting 
/ = 1 we recover the "harmonic" slicing condition, while, 
by setting / = q/a, with q an integer, we recover the gener- 
alized "1+log" slicing condition |31]. In particular, all the 
simulations discussed in this paper are done using condition 
(12.10b with / = 2 /a. This choice has been made mostly be- 
cause of its computational efficiency, but we are aware that 
"gau ge p athologies" could develop with the "1+log" slic- 
ings i 32l[33ll . 

As for the spatial gauge, we use one of the "Gamma-driver" 
shift conditions proposed in ll34ll . that essentially acts so as 
to drive the F* to be constant. In this respect, the "Gamma- 
driver" shift conditions are similar to the "Gamma-freezing" 
condition df T'^ — 0, which, in turn, is closely related to the 
well-known minimal distortion shift condition jisll . 

In particular, all the results reported here have been ob- 
tained using the hyperbolic Gamma-driver condition, 

d^f3' =Fdtr-r^dtf3\ (2.11) 

where F and r] are, in general, positive functions of space 
and time. For the hyperbolic Gamma-driver conditions it is 
crucial to add a dissipation term with coefficient 77 to avoid 
strong oscillations in the shift. Experience has shown that by 
tuning the value of this dissipation coefficient it is possible to 
almost freeze the evolution of the system at late times. We 
typically choose F — ja and ij — 2 and do not vary them in 
time. 



B. Evolution of the hydrodynamics equations 

The stellar models are here treated in terms of a perfect fluid 
with stress-energy tensor 

Tt"" = phu'^u" ^pg^"" , (2.12) 

h = ! + €+?- , (2.13) 
P 

where h is the specific enthalpy, e the specific internal energy 
and p the rest-mass density, so that e = p{l + e) is the en- 
ergy density in the rest-frame of the fluid. The equations of 
relativistic hydrodynamics are then given by the conservation 
laws for the energy, momentum and baryon number 

VuT'^" = , 

(2 14) 

V^ipu^^) = , ' ■ ^ 

once supplemented with an EOS of type p = p{p, e). While 
the code has been written to use any EOS, all the simulations 
presented here have been performed using either an isentropic 
"polytropic" EOS 

p^Kp^, (2.15) 

where K is the polytropic constant and F the adiabatic expo- 
nent, or a non-isentropic "ideal-fluid" (F-law) EOS 

p=(r-l)pe. (2.16) 

Note that, with the exception of the polytropic EOS ( 12.151 1. the 
entropy is not constant and thus the evolution equation for e 
needs to be solved. 

An important feature of the Whisky code is the imple- 
mentation of a conservative formulation of the hydrodynamics 
equations in which the set of equations ( 12.141 ) is written in a 
hyperbolic, first-order and flux-conservative form of the type 

atq + 9,fW(q) = s(q), (2.17) 

where f (q) and s(q) are the flux-vectors and source terms, 
respectively (see ref. |36] for an explicit form of the equa- 
tions). Note that the right-hand side (the source terms) de- 
pends only on the metric, and its first derivatives, and on the 
stress-energy tensor In order to write system (12.14b in the 
form of system ( 12.171 ). the primitive hydrodynamical variables 
(i.e. the rest-mass density p and the pressure p (measured in 
the rest-frame of the fluid), the fluid three-velocity (mea- 
sured by a local zero-angular-momentum observer), the spe- 
cific internal energy e and the Lorentz factor W ~ au'^) are 
mapped to the so called conserved variables q = {D,S'^,t) 
via the relations 

D = ^Wp, 

S' = ^phW^v' , (2.18) 

r = V7 {phW^ -p)-D . 

As previously noted, in the case of the polytropic EOS ( 12.151 ). 
one of the evolution equations (namely the one for r) does not 
need to be solved as the internal energy density can be readily 
computed by inverting the relation ( 12.161 ). Additional details 
of the formulation we use for the hydrodynamics equations 
can be found in | 36i] . 
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FIG. 1: Position on the {M/R,j3) plane of the stellar models con- 
sidered. Indicated respectively with stars and filled circles are the 
stable and unstable models belonging to a sequence of constant rest 
mass Mo — 1.51 Mq. Open circles refer instead to models which 
do not belong to the sequence, are also unstable and were first inves- 
tigated in ref. IT^ . Finally, the inset shows a magnification of the 
region where the threshold of the instability has been located for the 
sequence of models investigated. 



III. INITIAL DATA 

The initial data for our simulations are computed as sta- 
tionary equilibrium solutions for axisymmetric and rapidly 
rotating relativistic stars in polar coordinates 1371] , In gen- 
erating these equilibrium models we assumed that the metric 
describing an axisymmetric and stationary relativistic star has 
the form 



^ei'+'^dt^ + e^-^r^ sin^ e{d(j) - udtf 



+e^^{dr^ +r^de^) 



(3.1) 



where /i, v, uj and ^ are space-dependent metric functions. 
Similarly, we assumed the matter to be characterized by a non- 
uniform angular velocity distribution of the form 



Ik 



(il - w)r2 sin^ 



'1V 



1 - (r2 



2 a^-iv 



(3.2) 



where r^, is the coordinate equatorial stellar radius and the co- 
efficient A is a measure of the degree of differential rotation, 
which we set to A = 1 in analogy with other works in the lit- 
erature. Once imported onto the Cartesian grid and through- 
out the evolution, we compute the angular velocity 51 (and the 
period P) on the (x, y) plane as 



^ cos ( 
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(3.3) 



and other characteristic quantities of the system, such as the 
haryonic mass il/g, the gravitational mass A/, the angular mo- 
mentum J, the rotational kinetic energy T and the gravita- 
tional binding energy W as 



M 
Mo 

J 

T 
W 



J d\ (-2T° + T;) aV7 , (3.4) 

(3.5) 



d-'xD 



d\ De , 



\ Jd^xnr^a^ , 
T + Eint + Mo - M , 



(3.6) 

(3.7) 

(3.8) 
(3.9) 



where a y/j is the square root of the four-dimensional metric 
determinant. We recall that the definitions of quantities such 
as J, T, W and f3 are meaningful only in the case of stationary 
axisymmetric configurations and should therefore be treated 
with care once the rotational symmetry is lost. 

All the equilibrium models considered here have been cal- 
culated using the relativistic polytropic EOS ( |2.15t with K = 
100 and F = 2 and are members of a sequence having a con- 
stant amount of differential rotation with A ~ 1 and a constant 
rest mass of A/q ~ 1.51 A/q (a part of this sequence has also 
been considered in refs. llssll and ll39ll as models A8 — AlO). 
These are collected in Fig. [T| in a compactness/instability- 
parameter plot, where we have indicated with stars the stable 
models (S1-S8) and with filled circles the unstable ones (Ul- 
Ul 1). In addition, in order to compare with previous results, 
we have also considered three other models (D2, D3, D7), 
first investigated in ref. [T9^, which have larger masses and 
compactnesses. These models are unstable and are marked by 
circles in Fig. [T| Finally, the inset shows a magnification of 
the region where the threshold (indicated with a dashed line) 
of the instability has been located for the members of the se- 
quence. 

The main properties of all the considered models are re- 
ported in Table |I] while we show in Fig. |2] the profiles of 
the rest-mass density p (left panel) and of the rotational an- 
gular velocity 51 (right panel) for some of the models in the 
constant-rest-mass sequence. Note that the position of the 
maximum of the rest-mass density coincides with the center 
of the star only for models with low (3; for those with a larger 
/3, the maximum of the rest-mass density resides, instead, on 
a circle on the equatorial plane. Finally, indicated with a dot- 
dashed line in Fig.|2]is the profile for the first unstable model 
(Ul) with (3 = 0.255. Note that this is not the first model 
having an off-centered maximum of the rest-mass density. 

As mentioned in the Introduction, numerical simulations 
of the dynamical bar-mode instability have traditionally been 
sped up by introducing sometimes very large initial m=2 de- 
formations. The rationale behind this is simple since a seed 
perturbation has the effect of reducing the time needed for 
the instability to develop and thus the computational costs. 
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FIG. 2: Initial profiles of the rest-mass density p (left panel) and of the angular velocity Q. (right panel) for models S8, S7, S2, Ul, U3, Ull 
and U13. Indicated with a dot-dashed line is the profile for the first unstable model (Ul) with 13 = 0.255. Note that this is not the first model 
having an off-centered maximum of the rest-mass density. 



However, as we will discuss in detail in Section IVl CI the in- 
troduction of any perturbation (especially when this is not a 
small one) may lead to spurious effects and erroneous inter- 
pretations. Although in almost all of our simulation we have 
evolved purely equilibrium models and simply used the trun- 
cation errors to trigger the instability, we have also considered 
models which are initially perturbed so as to determine the ef- 
fect of these perturbations on the evolution of the instability. 
In these cases, we have modified the equilibrium rest-mass 
density po with a perturbation of the type 



rV. METHODOLOGY OF THE ANALYSIS 

A number of different quantities are calculated during the 
evolution to monitor the dynamics of the instability. Among 
them is the quadrupole moment of the matter distribution, 
which we compute in terms of the conserved density D rather 
than of the rest-mass density p or of the Too component of the 
stress energy momentum tensor 



(4.1) 



Sp2ix,y,z) = 62PQ- 



(3.10) 



where S2 is the magnitude of the m=2 perturbation (which 
we usually set to be 62 — 0.01 — 0.3). This perturbation has 
then the effect of superimposing on the axiaUy symmetric ini- 
tial model a bar-mode deformation that is much larger than 
the (unavoidable) m=4-mode perturbation introduced by the 
Cartesian grid discretization. In addition to a bar-mode defor- 
mation and in order to test the effect of a pre-existing m=l- 
mode perturbation we also used in some cases {cf. Section 
IVlCb an TO=l-mode density perturbation of the type 



Spi{x,y,z) =Sisin[(f> + 



2Trw 



Po 



(3.11) 



with Si = 0.01. Finally, after the addition of a perturbation of 
the type ( 13. 10b or ( 13.111 ), we have re-solved the Hamiltonian 
and momentum constraint equations, in order to enforce that 
the initial constraint violation is at the truncation-error level. 



Of course, the use of D in place of p or of Too is arbitrary 
and all the three expressions have the same Newtonian limit. 
However, we prefer the form ( 14.11 ) because D is a quantity 
whose conservation is guaranteed by the form chosen for the 
hydrodynamics equations ( 12.171 ). The time variation of ( 14.11 ) 
(or, rather, suitable combinations of its second time deriva- 
tives) will then be used in Section IVlllI to characterize the 
gravitational-wave emission from the instability. 

Once the quadrupole moment distribution is known, the 
presence of a bar and its size ma y be usefully quantified in 
terms of the distortion parameters 1I2III 



jxx _ jyy 
Jxx _j_ Jyy 

2 py 

Jxx _j_ Jyy 



(4.2) 
(4.3) 
(4.4) 



In addition, the quantity (14. 2t can be conveniently used to 
quantify both the growth time of the instability and the 
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Model 



U13 
U12 
Ull 
UIO 
U9 
U8 
U7 
U6 
U5 
U4 
U3 
U2 
Ul 



SI 
S2 
S3 
S4 
S5 
S6 
S7 
S8 



D2 
D3 
D7 



(10-4) 



0.5990 
0.9940 
1.0920 
1.1960 
1.2840 
1.3470 
1.4060 
1.4810 
1.5530 
1.5880 
1.6720 
1.7230 
1.8120 



0.20010 
0.24150 
0.25010 
0.25860 
0.26550 
0.27030 
0.27470 
0.28030 
0.28560 
0.28810 
0.29430 
0.29780 
0.30500 



1.8600 
1.8850 
1.9160 
1.9620 
2.1280 
2.2610 
2.7540 
3.8150 



0.30700 
0.30900 
0.31100 
0.31500 
0.32600 
0.33600 
0.37040 
0.44370 



3.1540 
3.7250 
2.7960 



0.27500 
0.30000 
0.30000 



Re 



Mo 



M M/Re 



24.31 
23.52 
23.31 
23.08 
22.88 
22.73 
22.59 
22.40 
22.22 
22.13 
21.92 
21.78 
21.54 



1.505 
1.508 
1.507 
1.508 
1.508 
1.508 
1.509 
1.508 
1.508 
1.508 
1.506 
1.508 
1.499 



1.462 
1.462 
1.460 
1.460 
1.460 
1.460 
1.460 
1.459 
1.458 
1.458 
1.456 
1.457 
1.448 



0.0601 
0.0622 
0.0627 
0.0633 
0.0638 
0.0642 
0.0647 
0.0651 
0.0656 
0.0659 
0.0664 
0.0669 
0.0672 



21.42 
21.35 
21.27 
21.14 
20.70 
20.32 
19.03 
16.70 



1.512 
1.510 
1.512 
1.504 
1.510 
1.505 
1.506 
1.506 



1.460 
1.458 
1.459 
1.452 
1.456 
1.449 
1.447 
1.439 



0.0682 
0.0683 
0.0686 
0.0687 
0.0703 
0.0713 
0.0760 
0.0862 



18.29 
17.85 
19.56 



2.771 
2.640 
2.188 



2.587 
2.466 
2.075 



0.1414 
0.1382 
0.1061 



Pa 

(ms) 



(ms) 



3.747 
3.591 
3.541 
3.496 
3.457 
3.428 
3.402 
3.363 
3.326 
3.310 
3.261 
3.241 
3.164 



1.723 
1.599 
1.572 
1.542 
1.517 
1.500 
1.484 
1.465 
1.446 
1.437 
1.417 
1.404 
1.386 



3.910 
3.654 
3.597 
3.536 
3.486 
3.450 
3.417 
3.377 
3.339 
3.321 
3.279 
3.251 
3.214 



3.191 
3.170 
3.160 
3.111 
3.050 
2.965 
2.741 
2.322 



1.368 
1.364 
1.356 
1.348 
1.308 
1.282 
1.189 
1.048 



3.180 
3.170 
3.153 
3.137 
3.057 
3.002 
2.812 
2.531 



7.620 
6.827 
5.386 



0.735 
0.731 
0.959 



2.051 
2.026 
2.442 



T 



W 



2.183 
2.272 
2.284 
2.302 
2.316 
2.325 
2.334 
2.341 
2.346 
2.351 
2.352 
2.360 
2.336 



7.764 
8.228 
8.327 
8.461 
8.575 
8.659 
8.741 
8.832 
8.920 
8.970 
9.061 
9.146 
9.167 



2.384 
2.378 
2.385 
2.363 
2.386 
2.369 
2.360 
2.255 



9.388 
9.396 
9.458 
9.439 
9.736 
9.859 
10.56 
11.96 



9.256 
8.450 
5.390 



35.28 
33.11 
21.06 



0.2812 
0.2761 
0.2743 
0.2721 
0.2701 
0.2686 
0.2671 
0.2651 
0.2631 
0.2621 
0.2596 
0.2581 
0.2549 



0.2540 
0.2531 
0.2522 
0.2503 
0.2451 
0.2403 
0.2234 
0.1886 



0.2624 
0.2544 
0.2561 



TABLE I: Main properties of the stellar models used in the simulations. Starting from the left the different columns report: the central rest- 
mass density pc, the ratio between the polar and the equatorial coordinate radii r-p/re, the proper equatorial radius Re, the rest mass Mq, the 
gravitational mass M, the compactness M / Re, the total angular momentum J, the rotational periods at the axis Pa and at the equator Pe, the 
rotational energy T and the binding energy W , and their ratio 13 (instability parameter). 



oscillation frequency of the unstable bar once the insta- 
bility is fully developed. In practice, we perform a nonlinear 
least-square fit of the computed distortion 77+ {t) with the trial 
function 

77+ (t) = e*/^B cos(2^ /3 t + (/-o) . (4.5) 

Note that all quantities ( I4.2b - (l4.4b are expressed in terms of 
the coordinate time t and do not represent therefore invariant 
measurements at spatial infinity. However, for the simulations 
reported here, the lengthscale of variation of the lapse function 
at any given time is always larger than twice the stellar radius 
at that time, ensuring that the events on the same timeslice are 
also close in proper time. 

Unless stated differently, we generally do not impose any 
boundary condition enforcing certain symmetries. As a result, 
during the evolution the compact star is not constrained to be 
centered at the origin of the coordinate system and, in order 
to monitor the relative motion of the rest-mass density distri- 
bution with respect to the coordinate system, we compute the 
first momentum of the rest-mass density distribution 

Xl^ = L (d\px\ (4.6) 
M J 



where M = J (fix p. These quantities are reminiscent of the 
Newtonian definition of the centre of mass of the star but, be- 
cause they are not gauge-invariant quantities, they are not ex- 
pected to be constant during the evolution. However, since in 
a Newtonian framework a time-variation of one of the 
would signal a nonzero momentum in that direction, we mon- 
itor these quantities as a measure of the overall accuracy of the 
simulations. Note also that, since the concept of the centre of 
mass is well defined in a Newtonian context only, equivalent 
definitions to ( I4.6l l could be made in terms of D or of Tqq. 
We have verified that in our simulations no significant quan- 
titative differences are present among the possible alternative 
definitions. 

In addition, as a fundamental tool to describe and under- 
stand the nonlinear properties of the development and satu- 
ration of the instability, we decompose the rest-mass density 
into its Fourier modes so that the "power" of the m-ih mode 
is defined as 



P„, = dhpe'"^^ (4.7) 
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FIG. 3: Mode-phases (solid line for the m=2 mode, dashed for the 
m=4 mode and dot-dashed for the m=\ mode) at different values of 
w overlapped with isocontours of the rest-mass density for model 
Ull at 25.7ms. 

and the "phase" of the m-th mode is defined as 

= arg(P™) . (4.8) 

The phase <?!>m essentially provides the instantaneous orienta- 
tion of the m-th mode when this has a nonzero power and is 
expected to have a harmonic time dependence when the cor- 
responding mode has a fully developed mode-component. 

An important clarification to make is that, despite their de- 
nomination, the Fourier modes ( 14. 7t do not represent proper 
eigenmodes of oscillation of the star. While, in fact, the latter 
are well defined only within a perturbative regime, the former 
simply represent a tool to quantify, within a fully nonlinear 
regime, what are the main components of the rest-mass distri- 
bution. Stated differently, we do not expect that quasi-normal 
modes of oscillations are present but in the initial and final 
stages of the instability, for which a perturbative description 
is adequate. 

Note also that the diagnostic quantities (14.71 ) are closely re- 
lated to the "dipole diagnostic " D = Pi /A 'l and "quadrupole 
diagnostic" Q = P2/M of ref. 1,17, 1. For some selected mod- 
els we have restricted the integration domain in eqs. ( 14.71 ) and 
( I4.8I 1 to the equatorial [i.e. {x, y)] plane and performed an 
integration in the azimuthal angle (f> only. In this way the cor- 
responding quantities 

P,n{nj) = I d(?!> p(tucos(0),Ti7sin(0))e™'^ , (4.9) 

0„(tu) = arg{P^{njj) , (4.10) 

have an explicit dependence on the cylindrical radial coordi- 
nate U7 only. The quantities ( 14.101 ) have the advantage that 
they can be used to check the "coherence of the mode" since 



0„j (tu) should be independent of tu when the m-th mode is a 
global property of the matter distribution. As an example we 
show in Fig.[3]the phases for the m=l, 2 and m=4 modes for 
model Ull when the bar is still fully developed, just before 
the bar loses its coherence. Note that the m=\ mode shows a 
spiral-like pattern inside the star, while both ^2 and 04 acquire 
a radial dependence in the outer parts of the star, where the bar 
deformation is absent. A similar behaviour for the (t>m (07) has 
been observed in all the performed simulations. 

V. GENERAL FEATURES OF THE DYNAMICS 
A. The tests on stable models 

Before investigating the nonlinear dynamics of unstable 
stellar models, we have caiTied out a systematic investigation 
of the ability of our code to perform long-term stable and ac- 
curate evolutions of stable stellar models. In particular, we 
have considered the time evolution of two of the differentially 
rotating models discussed in ref. ll38l[39tl . namely models S7 
and S8, and have followed their dynamics for 24 and 35 axial 
rotation periods, respectively. In both cases the stellar mod- 
els remain stable and the density and velocity fluctuations in 
the stellar interior are smaller than 2% during the whole sim- 
ulation. This is a rather remarkable result in fully 3D simula- 
tions and it is worth stressing that the simulations reported in 
ref. 1 40] were not able to go beyond 3 orbital periods for sim- 
ilar values of the grid size and spacing (we recall that in [4(J 
a second-order TVD method with the MC limiter was used in 
place of the third-order PPM method used here). 

In addition, for a more quantitative check of the accu- 
racy of our simulations, we have computed the frequency of 
the /-mode using the normalized power spectrum (Lomb's 
method (a^) of the coordinate time evolution of the central 
rest-mass density. The calculated values of 791 Hz for model 
S8 (A9) and of 674 Hz for model S7 (AlO) are in very good 
agreement, with the values of 809 Hz and 685 Hz reported 
in ref. [39J and computed using a 2D grid in spherical coor- 
dinates but in the conformally flat approximation of General 
Relativity. 

B. Common feature of unstable models 

In this subsection, we discuss some of the general features 
of the dynamics of unstable models, postponing to the fol- 
lowing sections the discussion of more detailed aspects of the 
instability. Here we will focus in particular on the dynam- 
ics of three representative unstable models, namely U3, Ull 
and U13, which have been selected so that their increasing 
values for the (3 parameter cover the whole range of inter- 
est. For these simulations, we have used a spatial resolution 
Ax = 0.5 Mq and a grid of 193 x 193 x 68 cefls and imposed 
a reflection symmetry with respect the {x, y) plane. As a re- 
sult, between 80 and 90 gridpoints cover the stars along the x 
and y axes at time i = 0. Note that all the simulations reported 
here make use of a uniform grid with the location of the outer 
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FIG. 4: Snapshots of the evolution of models U3, Ull and U13 at various times. The different columns refer to the three models and show 
isodensity contours for p = 0.9, 0.8, 0.7, 0.6, 0.5^ -' X Pmax> where j — 1, . . . , 6 and pmax is the maximum value of p in each panel. The 
above models were evolved on a 193 x 193 x 68 grid with grid coordinate resolution of 0.5 Mq (0.74km) and imposing equatorial symmetry. 
The time evolution of some quantities characterising these models is reported in Figs.|5][6]and|7] 
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FIG. 5: Time evolution of the instability for model U3. The top panel 
shows the behaviour of the quadrupole distortion parameter 77+ [c/ 
eq. J4.4b 1. the middle panel reports the behaviour of the power in the 
Fourier modes m=l, 2, 3 and 4, while the bottom panel displays the 
phase of the m=2 mode. 



boundary being rather close to the stellar surface; this makes 
the extraction of gravitational waves difficult and accounts for 
a very small but nonzero loss of mass and angular momentum 
(because of matter escaping the computational box). 

In Fig.|4]we show some representative snapshots of the rest- 
mass density at four different times for three different unstable 
models (one column for each model). In particular, each row 
refers to one of the four representative stages in which the dy- 
namics of the bar can be divided. These are: {a) exponential 
growth of the to=2 mode and m=3 mode (first row); {h) sat- 
uration of the instability, development of spiral arms and pro- 
gressive attenuation of the bar deformation (second row); (c) 
crossing of to=3 mode and m=4 mode and consequent attenu- 
ation of the bar, emergence of the m=l mode as the dominant 
one (third row); {d) suppression of the bar deformation and 
emergence of an almost axisymmetric configuration (fourth 
row). Note that while these stages are present in all these 
three models, the coordinate times at which they take place 
(indicated in the upper part of each panel), as well as the am- 
plitude of the deformation, depend on the parameters defining 
the initial models, most notably /3 and Mq. 

Understanding the occurrence of these four stages during 
the onset, development and suppression of the bar deforma- 
tion represents our effort to go beyond the standard phe- 
nomenological discussion of the nonlinear dynamics of the 
instability often encountered in the literature. An important 
tool in this discussion will be offered by the time evolution of 
the Fourier mode-decomposition (|4.7| i discussed in Sec. |IV] 
As we will show below, relating the evolution of these quanti- 
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FIG. 6: The same as Fig.[5]but for model Ull. 



ties to the evolution of the mode phases and to the changes 
in the deformation of the star »7+, 77x will allow us to provide 
a consistent description of the four stages of the instability. 

We start our discussion by reporting in Figs. |5] |6] and |7] 
the history of the instability for models U3, Ull and U13. 
Starting from the upper panels, these figures show: the time 
evolution of the distortion parameter (a very similar be- 
haviour can be shown for the other distortion parameter Tyx ), 
the power P„i in the first four m-modes and the evolution of 
the phase of the m=2 mode. Note that at the beginning of 
the simulation, as a result of the Cartesian discretization, the 
TO=4 mode has the largest power. While this can be reduced 
by increasing the resolution, the m=4 deformation plays no 
major role in the development of the instability, which is soon 
dominated by the lower-order modes. 

The initial phase of the instability [stage fa) in the previ- 
ous classification] is clearly characterized by the exponen- 
tial growth of the m=2 mode and m=3 mode, the latter one 
having a smaller growth rate. A first interesting mode cou- 
pling takes place when the exponentially growing m=3 mode 
reaches the same power amplitude of the m=4 mode, at which 
point the two modes exchange their dynamics, with the m=4 
mode growing exponentially and the m=3 mode reaching sat- 
uration. At approximately the same time, the m=l mode also 
starts to grow exponentially but with a growth rate which is 
smaller than that of the other modes. Note that this "mode- 
amplitude crossing" between the to=3 and m=4 modes also 
signals the time when collective phenomena start to be fully 
visible (note that this mode-amplitude crossing is distinct 
from the "avoided-crossing" observed when studying mode 
eigenfrequencies along sequences of stellar models). This is 
shown with the first vertical dotted line in Figs. |5] |6]and|7] 
marking the time when the distortion parameter starts being 
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FIG. 7: The same as Fig.|5]butformodelU13. 



FIG. 8: Schematic evolution of the collective modes [Eg. 14.71 of the 
rest-mass density p. In this diagram the instability is distinguished 
in four representative stages: (a) exponential growth of the m=2 
mode and m=3 mode; (b) saturation of the instability, development 
of spiral arms and progressive attenuation of the bar deformation; 

(c) crossing of m=3 mode and m=4 mode and consequent attenua- 
tion of the bar, emergence of the m=l mode as the dominant one; 

(d) suppression of the bar deformation and emergence of an almost 
axisymmetric configuration. 



appreciably different from zero (upper panel) and the m=2 
phase assumes a harmonic time dependence (lower panel). 
This stage continues until the m=2 mode reaches its maximum 
power and the bar has reached its largest extension. During the 
following phase [stage (b)] the bar instability has reached a 
nonlinear saturation, accompanied by the development of spi- 
ral arms which are responsible for ejecting a small amount of 
matter and for a progressive attenuation of the bar extension 
(see discussion in Sec. IVl Ab . Furthermore, when the expo- 
nentially growing m=\ mode reaches the same power ampli- 
tude of the m=3 mode, the latter, whose growth had slowed 
down for a while, returns to grow exponentially. 

The following phase of the instability [stage fcj] sees modes 
m=l, 3 and 4 reach comparable powers and this marks the 
time when the bar deformation has a sudden decrease. As 
a result of this crossing among the three modes, only the 
m=\ mode will continue to grow, while the to=3 and the 
m=A modes are progressively damped. Finally, stage (d) 
starts when the growing m=l mode reaches power amplitudes 
comparable with those of the to=4 mode and the final mode- 
amplitude crossing takes place. This marks a distinct loss of 
the bar deformation and the emergence of an almost axisym- 
metric rapidly rotating star. This is shown with the second ver- 
tical dotted line in Figs. |6] and |7] highlighting when the distor- 
tion parameter is significantly reduced (upper panel) and the 
TO=2-mode phase loses its harmonic time dependence (lower 
panel). A schematic and qualitative diagram summarizing the 
evolution of the power in the first four m-modes as discussed 
above is shown in Fig. [8] and can be used as an aid for the 
interpretation of the quantities computed in Figs.|5] |6]and|2l 

We note that the lack of a perturbative study of this process 
beyond the linear regime leaves the origins of this interaction 
between modes still unclear. Furthermore, since the growth 



of the m=l mode is not clearly exponential, especially for 
slightly overcritical models {cf. Fig.|5]for model U3), we have 
referred to this process as "mode coupling" rather than con- 
sidering it as the evidence of an m=l instability. Additional 
perturbative work in this respect will help clarify this aspect. 

The general and common features of the dynamics of the 
bar-mode instability as deduced from the numerical simula- 
tions can be summarized as follows: 

• the bar deformation is, in general, not a persistent phe- 
nomenon but is suppressed rather rapidly and over a 
timescale which is of the order of the dynamical one 
(see also the following Section for an additional discus- 
sion on this); 

• nonlinear mode couplings take place during the evolu- 
tion and these allow for the growth of other modes be- 
sides the fastest-growing m=2 mode; 

• the growth of other modes has the overall impact of 
progressively attenuating the m=2 mode and, conse- 
quently, the bar deformation, after the instability has 
saturated; 

• for slightly supercritical models {e.g. U3), when the 
power amplitude of the m=\ mode has become com- 
parable with the one in the m=2 mode, the bar defor- 
mation is suppressed and the star evolves towards an 
almost axisymmetric configuration; 

• for largely supercritical models {e.g. U13), the dynam- 
ics of the instability are so violent and the stellar model 
so far from equilibrium that the strong bar deformation 
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is lost even in the absence of mode-coupling effects (see 
discussion in Sec lVIBl l. 



VI. DETAILED FEATURES OF THE DYNAMICS 

In this Section we discuss some detailed aspects of the in- 
stability, concentrating our attention on the impact that differ- 
ent values of f3, different boundary conditions, different val- 
ues of the initial perturbations, different EOSs and different 
grid resolutions or boundary locations have on the onset and 
development of the instability. We note that while many of 
these different prescriptions do not induce qualitative changes, 
some of them do change the initial relative amplitude of the 
different modes [and hence the simulation time needed for 
the instability to develop and the orientation of the bar in the 
{x,y) plane at a given time during the instability]. In these 
cases, in order to make meaningful comparisons, we remove 
these offsets by choosing a suitable shift in time At and in 
phase Aip in such a way that the distortion parameters of the 
reference model ij!^ and of the new one 77+ have the maximal 
overlap and are related as 

r^^fit) ~ ar^+{t + At) + p-q^ {t + At) , (6.1) 
where a — cos(A(/)), /3 = sin(A0). 

A. Dependence on f3 

The parameter f3 plays a very important role in determining 
the properties of the nonlinear dynamics of the instability both 
with regard to the growth rate T3 and to the duration t-^ of 
the saturation stage [stage (b) of Fig. [8). While the relation 
between (3 and will be discussed in more detail in Sec. lVIII 
we here concentrate on how the dynamics of the bar, once 
formed, depend on the degree of overcriticality, i.e. on /3 — /3c, 
where (3c marks the separation between stable and unstable 
models. To illustrate this we will consider two models which 
are representative of the whole set considered in Table |T] and 
which have very different values of /? and consequently very 
distinct behaviours: the largely overcritical model U13 and 
the slightly overcritical model U3. 

Model U13 has (3 = 0.2812 and is the most unstable of the 
studied models, since equilibrium models with larger (3 can- 
not be produced along the constant-rest-mass sequence cho- 
sen here. The clearest feature shown in Fig. |7]for this model 
is its very rapid growth rate (almost three times larger than the 
one for U3 as reported in Tab. but also its very effective 
suppression of the bar deformation. Indeed, the evolution is 
so rapid that it is very little affected either by initial perturba- 
tions or by the imposition of additional symmetries (see next 
subsections). In this case, in fact, the crossing between the 
m=2 mode and m=l mode is less evident and the bar defor- 
mation goes through large variations, as shown by the large 
oscillations in P2 after t ~ 16 ms, during which time the bar 
seems to disappear and then form again soon after. At about 
i ~ 20 ms the bar deformation starts disappearing in coinci- 
dence with the mode-amplitude crossing. Again as a result of 
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FIG. 9: The role of the 7r-symmetry on the dynamics for model U3. 
Shown from the top are the deformation parameter 77, the power in 
the m=2 mode and in the m=l mode (dashed line), the evolution of 
the position of the "centre of mass" (the horizontal dashed lines mark 
the edges of the central cell) and that of the rest mass. The dotted and 
continuous lines refer to simulations without and with 7r-symmetry, 
respectively. 



the very violent dynamics, the saturation stage is rather short 
and the bar is essentially lost after about 8 ms. 

Model U3, on the other hand, has /? = 0.2596 and shows 
dynamics which are in many respects the opposite of the ones 
discussed for model U13. As shown in Fig.|5] the mode evolu- 
tion is very smooth and, once formed, the bar persists without 
significant losses in power. The growth rate is clearly smaller 
and the stage of saturation is much longer (about 30 ms) and 
the growth of the m= 1 mode plays a major role in the damping 
of the 171=2 mode. The transition that leads to the disappear- 
ance of the bar is smooth and it requires many rotation peri- 
ods. Differently from model U13, in this case, the properties 
of the bar dynamics in the first stage are sensitive to the use 
of perturbations or to the imposition of additional symmetries 
(see next subsections). 

Overall, it is reasonable to expect that the persistence of the 
bar is strictly related to the degree of overcriticality, with the 
duration of the saturation t-^ tending to the radiation-reaction 
timescale for a model with (3 — (3c and to zero for a model 
with /? ^ (3c, for which the excess of kinetic rotational en- 
ergy may well produce a rapid disruption of the star (see Fig. 
fTTT i. The numerical values for have been estimated through 
a nonlinear fit to the evolution of P2 with 3 separated single 
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FIG. 10: The same as in Fig.|9]but for model U13. 



exponential functions in the three intervals [2, to], [ta, h] and 
[tb, tc], where ta and tb are two free parameters and tc marks 
the end of the simulation. The estimates for reported in 
Table are still too sparse to be able to delineate its depen- 
dence, beyond the evidence that oc |/3 — Pc\^'", where n is 
a positive number. Furthermore, the reported values have an 
error of about 1 ms, as a result of the used fitting procedure. 



B. The role of symmetries 

As mentioned in the Introduction, the issue of the persis- 
tence of the bar deformation has been rather controversial over 
the years. 

While previous calculations carried out in Newtonian 
physics and in the absence of symmetries have highlighted 
that the bar deformation can be rapidly suppressed as a re- 
sult of the growth of an rn=l-mode deformation lUt], subse- 
quent studies have attributed the growth of the odd mode to 
inaccurate numerical methods and supported the idea that the 
bar should be persistent over a radiation-reaction timescale 
and that the use of suitable symmetry conditions that remove 
the growth of the odd mode provides a more realistic de- 
scription of the bar dynamics [9]. In addition, it has been 
argued that once an 7Ti=2-mode perturbation has developed, 
only couplings with even modes should be expected and that 
the growth of any odd mode should therefore be considered a 
spurious numerical artifact. 
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FIG. II: Dynamics of the rotational kinetic energy T [cf. eq. l |3.8l )1 
and of the internal energy Eint [cf. eq. l |3.6H for models U3, UIl 
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We believe the above argument not to be valid, except in a 
linear regime and in the very idealized case in which it is pos- 
sible to inject exclusively an m=2 perturbation. In practice, 
however, any initial perturbation, either introduced ad hoc or 
by the truncation error, will excite both even and odd modes 
and all of these will couple once a nonlinear regime is reached. 

Having said this, it is nevertheless important to verify that 
the growth of the m=l mode detected in our simulations is not 
a numerical artifact (this is further discussed in Sec. VI.E) and 
that the argument made about the non-persistence of the bar 
deformation continues to hold also when boundary conditions 
with symmetries are introduced. For this reason, we evolved 
the models discussed in the previous Sections also with the 
use of the so-called 7r-symmetry, ensuring that /(nj, 0, z) = 
/(tn, 4> + TT, z) for any variable /(x*). Clearly, the presence of 
any odd mode is in this way impossible by construction. We 
report in Figs. |9lfT0l the results of the simulations for models 
U3 and U13 using this symmetry. The first two panels from 
the top show the deformation parameter rj{t) and the power in 
the 171=2 mode (solid line when the 7r-symmetry is enforced 
and dotted line otherwise) and in the m= \ mode (dashed line). 

Fig.|9]clearly shows that the bar deformation is essentially 
persistent in model U3 when the symmetry boundary condi- 
tions are applied, its power amplitude being just slowly atten- 
uated, mostly because of the entropy production via the non- 
isentropic EOS ( |2.16t and a possible small contribution due 
to the use of a tenuous atmosphere outside the star. However, 
it is important to note that, besides being convergent and sta- 
ble, the solution without 7r-symmetry is also accurate, as well 
as the one without 7r-symmetry. This is shown by the second 
panel from the bottom in Fig. |9l reporting the evolution of 
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FIG. 12: Effect of an initial m=2-mode perturbation on the dynamics of the deformation parameter ri{t) (top sub-panels) and of the modes 
P2{t) and Pi{t) (bottom sub-panels) for model U3 (left) and U13 (right), respectively. The continuous lines represent the evolution of the 
perturbed model after a suitable phase and time shifts. 



the position of the "centre of mass" as defined in eq. ( 14.6b . in 
the absence of 7r-symmetry. The two horizontal dashed lines 
in that panel mark the edges of the central cell and indicate 
that up to t ~ 50 ms the position of the centre of mass does 
not leave the central cell of the grid and that the exponential 
growth of the m=l mode (which becomes significant from 
t ~ 15 ms) cannot be related to a spurious numerical effect. 
After i ~ 50 ms the centre of mass starts to move away from 
the centre of the grid and also in this case the motion is not 
due to numerical accuracy but rather to the fact that a small 
amount of matter (about 2% of the initial one) is being lost 
from the grid as a result of the development of extended spiral 
arms. This is shown in the lower panel of Fig.|9] which reports 
the evolution of the rest mass when normalized to the initial 
value and which clearly shows that the motion of the centre 
of mass is related to the loss of rest mass through the grid and 
thus consistent with the conservation of linear momentum. As 
we will further discuss in Sec. lVIEl the mass loss and the con- 
sequent motion of the centre of mass can be reduced consider- 
ably when moving the outer boundary to larger positions. As 
mass leaves the grid, so does angular momentum, with losses 
that vary according to the model considered and ranging from 
~ 3% for model U3 up to ~ 20% for the more violent model 
U13. Note, however, that much smaller angular- momentum 
losses (i.e. less than 1% for all models) are in general mea- 
sured before the mass is shed across the computational bound- 
aries; as a result, angular momentum can be conserved to rea- 
sonable accuracy by using more distant outer boundaries (this 
improvement has been tried and a discussion on the changes 
introduced is presented in Sec. lVIEl ). 

Interestingly, the use of a 7r-symmetry does not produce 
a significant change in the case of model U13. This is true 



for the dynamics of the bar (cf. the solid and dotted lines in 
Fig. [TOb and also for the values of r^, and /g. We believe 
this is because the dynamics of this largely overcritical model 
are not dominated by the mode coupling, but rather by effi- 
cient conversion of rotational kinetic energy into internal en- 
ergy. As shown in Fig. [TT] which reports the time evolution 
of the internal and rotational kinetic energies when normal- 
ized to their initial values, model U13 experiences a dramatic 
and rapid increase in the internal energy at the expense of the 
kinetic one. This conversion of energy is the largest among 
the simulated models and so effective that mode-coupling ef- 
fects do not have time to develop. This explains why the use 
of symmetry conditions slightly reduces the attenuation of the 
bar, but cannot prevent its rapid disappearance. 

As a final remark we note that the use of a 7r-symmetry 
produces only small changes in the values of the bar-pattern 
frequency or of the growth time when compared with 
the corresponding values computed in the absence of symme- 
tries [cf. Table Hill. This is essentially because these boundary 
conditions do not alter the dynamics of the stage of exponen- 
tial growth of the bar. They can therefore be used to reduce 
the computational costs and pursue the systematic search for 
the threshold of the instability discussed in Sec. I VII 51 



C. The role of the initial perturbation 

Another common feature of previous works on the bar- 
mode instability has been the introduction of a sizeable m=2- 
mode perturbation of the type shown in e g. (13.101 1 with the 
goal of triggering the instabiUty S III [M M El III . This 
approach clearly reduces the computational costs but it is fully 
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FIG. 13: Effect of an initial perturbation on the dynamics of the deformation parameter rj (top panels) and of the mode powers P2{t), Pi{t) 
for model Ul 1. The left panel shows the effects of an initial m=2-mode perturbation, while the right one those of an m=l-mode perturbation. 



justified only when the triggered mode is the only unstable 
one. However, if other unstable modes exist, their develop- 
ment may be altered or even suppressed in the presence of a 
strong m=2-mode perturbation. This is particularly relevant 
for the analysis carried out in this work, which has pointed 
out that nonlinear mode couplings may trigger the growth of 
other modes and significantly modify the dynamics of the in- 
stability. 

We have therefore considered with care how the introduc- 
tion of an m=2-mode perturbation influences the onset and 
the development of the instability. More specifically, we 
have added an r7i=2-mode perturbation of the type shown in 
eq. (13.101 1 with 62 — 0.01, solved again the constraint equa- 
tions and observed that the impact this has on the develop- 
ment of the instability depends on the degree of overcriti- 
cality. In particular, for models near the threshold, such as 
U3, the perturbation does not induce changes in the saturation 
phase nor in the persistence of the bar, but it does have the 
effect of slightly altering the first part of the evolution, with 
an increase in the maximum distortion (which at saturation is 
^ 10% larger) as well as with an increase in the growth rate 
(the growth time Tg is reduced by ~ 10 %); see the left panel 
of Fig. [12] and Table Ull for a quantitative comparison. 

On the other hand, for models that are largely overcritical, 
such as U13, and in analogy with what discussed in the previ- 
ous Section for the use of symmetric boundary conditions, the 
introduction of a perturbation does not have a significant ef- 
fect and the dynamics are essentially unaltered (see the right 
panel of Fig. [T2] l. Finally, for models which are overcritical 
but not close to the threshold, such as Ull, the initial pertur- 
bation has a much smaller impact on both the growth rate and 
the maximum distortion (cf. Tablellll). but it does increase the 
duration of the bar This is due to the fact that the growth of 



the m=l mode is closely related to the one in the m=2 mode 
and its growth can be delayed and reduced if the latter has 
initially a non-negligible power Because of this, the time at 
which the two modes have comparable power will be different 
and in particular will be postponed in the perturbed case (see 
the left panel of Fig. \l3[ . Of course, the converse is also true 
and modified dynamics for this model are observed also when 
an m=l-mode perturbation of the type shown in eq. ( 13.10b is 
introduced with Si = 0.05. This is summarised in the right 
panel of Fig. [13] which shows that in this case the perturba- 
tion reduces the duration of the saturation stage. 

In summary, while the introduction of a seed perturbation 
(either in the form of an m=l mode or an m=2 mode) does 
not produce significant qualitative changes in the dynamics of 
the instability, it can result into quantitative changes, most no- 
tably in the growth rate, in the maximum distortion and in the 
persistence of the bar. The persistence of the bar, in particular, 
is enhanced when an TO=2-mode perturbation is present. The 
relevance of these results will need to be evaluated for those 
astrophysical scenarios in which long-lasting bars were simu- 
lated, but which were triggered through the introduction of a 
perturbation 1, 15. .23.1 . 



D. The role of the EOS 

Besides nonlinear mode coupling, another process that 
could in principle limit the persistence of the bar is the for- 
mation of shocks (either macroscopical or on smaller scales) 
that would convert the excess kinetic energy into internal one. 
In order to assess the importance of this process we have com- 
pared the evolution of the relevant unstable models when these 
are evolved using the non-isentropic EOS ( 12.16b and when us- 
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FIG. 14: Effects on the evolution of the deformation parameter 77 (top panels) and the modes P2{t) and P\{t) (bottom panels) of model Ull 
caused by: (left) the use of the adiabatic polytropic EOS; (right) the use of a larger simulation grid, where the distance of the outer boundary 
from the center is increased from 48 A/q to 66 Mq on the (x, 2/)-plane and from 32 Mq to 47 Mq along the z direction. 



ing the (isentropic) polytropic EOS ( |2.15t with K = 100 and 
r = 2. 

The results of this comparison are summarised for model 
Ull in Fig. [14] and indicate that the non-isentropic changes 
are indeed very small and that these do not produce any sig- 
nificant variations on the development of the instability and 
on the growth of the m=2 mode. Larger differences are seen 
in the growth of the m=l mode, but also these are very small 
and do not produce a qualitative change. Quantitative assess- 
ment of the changes produced by a different EOS are reported 
in Table [III but these are, overall, comparable with the error 
bar for the determination of and . Finally, all the consid- 
erations made here for model Ull apply also to models U13 
and U3, with model U3 being slightly more sensitive to the 
change in EOS {cf. Table Hill. These results indicate therefore 
that the effects of shock heating are likely to be unimportant at 
least for the development and evolution of the bar in isolated 
and old neutron stars. 



E. The role of grid spacing and size 

We finally report on the influence of the grid spacing and of 
the grid size on the development of the instability. To assess 
this we have performed several simulations of the interme- 
diate model Ull differing either in grid resolution or in the 
location of the outer boundaries. In particular, we have con- 
sidered grid resolutions Ax/Mq = 0.375, Ax/Mq — 0.5 
and Ax/Mq — 0.625 and found the code to be second-order 
convergent, with the coarsest resolution being just on the limit 
of the convergence regime (the results with Ax/Mq = 0.75 
are in fact not convergent at the expected rate). 



As summarized in Table we have found that the com- 
puted values of the instability parameters and do not 
vary significantly across the range of resolutions considered, 
with differences that are at most of about 3%. The same Table 
also contains information on the results obtained when com- 
paring simulations performed with Ax/Mq = 0.625 but with 
a larger computational domain, namely going from a compu- 
tational box, at this resolution, with extents [157 x 157 x 56] 
to one with extents [211 x 211 x 80]. Also in this case the 
changes in the dynamics are very small (see also the right 
panel of Fig. [T4b and essentially amount to a smaller loss of 
mass and angular momentum as some of the the matter in the 
spiral arms is thrown out of the computational grid (see also 
discussion in Sec lVIBl i. 



F. Comparison with previous studies 

To conclude this Section describing the dynamics of the in- 
stability, we comment on the important validation of the accu- 
racy of our simulations that comes from a comparison with re- 
sults previously published in the literature. We have focused, 
in particular, on the fully general-relativistic simulations pub- 
lished in ref. lfl9ll and repeated those relative to the models 
D2, D3 and D7 discussed there. These stars have instability 
parameters (3 rather close to the critical one, but are also more 
massive, with gravitational masses between 2 and 2.6 Mq, 
and have larger compactnesses {cf. TableUi. The development 
of the instability for one of these models (D2) is summarised 
in Fig. [15] and the computed distortion parameter is qualita- 
tively very similar to the one presented in ref. |19] and shows 
that in this compact star the bar [i.e. stage (b) of the classifi- 
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Model 


Ax/Mq 


<Si 


52 


EOS 


TT-sym 


grid size 


At 


ti 


t2 




/b 


V 


















(ms) 


(ms) (ms) 


(ms) 


(Hz) 


(max) 


(ms) 


U3 


0.500 


0.0 


0.0 


polytropic 


no 


medium 


0.76 


21 


26 


2.79 


552 


0.48 


13.8 


U3 


0.500 


0.0 


0.0 


ideal fluid 


no 


medium 





21 


26 


2.69 


547 


0.47 


12.9 


U3 


0.500 


0.0 


0.01 


ideal fluid 


no 


medium 


-15.94 


21 


26 


2.42 


548 


0.53 


13.0 


U3 


0.625 


0.0 


0.0 


ideal fluid 


yes 


medium 


-1.00 


21 


26 


2.82 


543 


0.43 


24.4 


U3 


0.625 


0.0 


0.01 


ideal fluid 


yes 


medium 


-16.28 


21 


26 


2.52 


547 


0.54 


24.5 


Ull 


0.500 


0.0 


0.0 


polytropic 


no 


medium 


5.35 


11 


14 


1.12 


497 


0.78 


8.6 


Ull 


0.500 


0.0 


0.0 


ideal fluid 


no 


medium 





11 


14 


1.15 


494 


0.78 


9.4 


Ull 


0.500 


0.0 


0.01 


ideal fluid 


no 


medium 


-8.55 


11 


14 


1.11 


494 


0.79 


9.9 


Ull 


0.375 


0.0 


0.0 


ideal fluid 


no 


medium 


1.64 


11 


14 


1.11 


492 


0.79 


10.5 


Ull 


0.625 


0.0 


0.0 


ideal fluid 


no 


medium 


1.79 


11 


14 


1.15 


492 


0.78 


9.0 


Ull 


0.625 


0.0 


0.0 


ideal fluid 


no 


large 


2.54 


11 


14 


1.15 


492 


0.78 


9.8 


Ull 


0.625 


0.0 


0.0 


ideal fluid 


yes 


medium 


1.39 


11 


14 


1.12 


494 


0.77 


13.8 


Ull 


0.750 


0.0 


0.0 


ideal fluid 


no 


medium 


3.79 


11 


14 


1.17 


493 


0.76 


10.8 


Ull 


0.625 


0.05 


0.0 


ideal fluid 


no 


medium 


6.20 


11 


14 


1.14 


495 


0.78 


6.6 


U13 


0.500 


0.0 


0.0 


polytropic 


no 


medium 


1.69 


10 


13 


0.94 


457 


0.86 


5.7 


U13 


0.500 


0.0 


0.0 


ideal fluid 


no 


medium 





10 


13 


0.95 


454 


0.85 


6.2 


U13 


0.500 


0.0 


0.01 


ideal fluid 


no 


medium 


-8.55 


10 


13 


0.93 


454 


0.86 


6.3 


U13 


0.625 


0.0 


0.0 


ideal fluid 


yes 


medium 


-0.16 


10 


13 


0.96 


453 


0.86 


6.5 


U13 


0.625 


0.0 


0.01 


ideal fluid 


yes 


medium 


-8.71 


10 


13 


0.94 


454 


0.86 


6.2 


D2 


0.500 


0.0 


0.0 


ideal fluid 


no 


medium 





9 


10.5 


0.90 


1053 


0.59 


- 


D2 


0.500 


0.0 


0.01 


ideal fluid 


no 


medium 


-6.54 


9 


10.5 


0.78 


1052 


0.67 




D2 


0.500 


0.0 


0.04 


ideal fluid 


no 


medium 


-7.57 


9 


10.5 


0.77 


1056 


0.67 




D3 


0.500 


0.0 


0.0 


ideal fluid 


no 


medium 





9 


12 


1.54 


1086 


0.38 




D7 


0.500 


0.0 


0.0 


ideal fluid 


no 


medium 





11 


14 


1.74 


821 


0.48 





TABLE II: Main properties of the initial part of the instability for the stellar models used in the simulations. Starting from the left the different 
columns report: the grid spacing Ax/Mq, the amplitudes of the initial perturbations in the m=l and m=2 modes 5i^2, the EOS, the symmetry 
and the grid size used, the time shift At [cf. eq. ( 16. IH , the times ti and t2 between which the growth-times and the frequencies are 
computed, the maximum value of the distortion parameter r/, and the duration of the bar deformation . 



cation made in Sec. lVBI is very rapidly attenuated as a result 
of the development of spiral arms, which are also responsible 
for a small loss of mass. 

The frequencies and the growth times found for these mod- 
els are also in good agreement with those reported in ref. lfl9ll . 
but are not identical; differences are of about 10% (see Ta- 
ble for a close comparison). While there are several dif- 
ferences in the numerical codes used, it should be noted that 
the simulations reported in ref. 1 19] made use of a substantial 
perturbation in the m=2 mode, with an equivalent 62 = 0.3. 
Although we were not able to reproduce exactly the dynam- 
ics of these models (no convergent solution of the constraint 
equations was found once such a large perturbation was intro- 
duced), we recall that large perturbations for models near the 
threshold do induce a change in the growth rates and effec- 
tively reduce the growth times (see discussion in Sec. lVICT l. 

We believe therefore that the use of a smaller or zero per- 
turbation is the largest source of the difference with the corre- 
sponding simulations in ref. iT9ll . which we can nevertheless 
reproduce to very good precision. 



VII. DETERMINATION OF THE THRESHOLD 

An important consequence of the high accuracy of our sim- 
ulations it that it has allowed for a rather precise determination 
of the stability threshold for the sequence of models consid- 
ered here. Of course, the determination of the threshold to the 
third significant figure has little but academic interest, as it is 
expected not to be universal, but rather to depend (although 
weakly) on properties such as the degree of differential rota- 
tion, the compactness, the EOS, etc.. Nevertheless, this is a 
useful exercise for at least two different reasons. Firstly, it 
helps in characterizing the dynamics of slightly supercritical 
models and, secondly, it allows for a direct comparison with 
perturbative studies, highlighting when the latter cease to be 
accurate and how they can be improved. 

In what follows we discuss two different methods we have 
used to this aim. The first one simply selects the initial model 
with the lowest value of /3 for which an exponential growth 
in the distortion parameter is seen. The second one, on the 
other hand, uses the classical Newtonian stability analysis of 
Maclaurin spheroids for incompressible self-gravitating fluids 
in equilibrium [IJ to extrapolate the position of the threshold 
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FIG. 15: The same as in Fig. [9] but for model D2. In addition, the 
bottom panel shows the evolution of the phase for the m =2 mode. 



also in a fully general-relativistic regime. Interestingly, the 
results of the two approaches agree to high precision. 



A. First method: dynamical evaluation 

The approach followed here is straightforward and consists 
in performing a number of simulations for models with de- 
creasing values of the instability parameter j3 and in determin- 
ing the critical value (3^ as the smallest one for which an expo- 
nential growth of the distortion parameter rj is observed. More 
specifically, we have performed a set of 10-ms simulations of 
models characterized by values of (3 between /3 = 0.24 and 
(3 — 0.255, interval in which the instability was reported to 
develop ITqI l20i, .2 1,1 . In order to remove a possible contam- 
ination by the Cartesian discretization through the power in 
the m=4 mode, all the models had a very small m=2-mode 
perturbation with 82 = 0.04, making this the largest mode 
power initially. Furthermore, since this kind of initial m=2- 
mode perturbation always generates, at least temporarily, a 
growth of the distortion, we classified as unstable those mod- 
els for which an exponentially growing bar deformation was 
observed for the whole simulated 10 ms. 



-0.1 



SI 

ui 




30 
t (ms) 

FIG. 16: Dynamics of the instability near the threshold. The two 
panels summarise two long-term evolutions for model SI (continu- 
ous line) and Ul (dotted line) and show the distortion parameter and 
the power in the m=2 and m=4 modes. 



As a result of this set of simulations and of additional re- 
finements of the bracketing interval for the instability, we have 
concluded that the threshold had to be found between models 
U2 and S2, although it was not yet obvious whether any or 
all of these models were unstable. A more precise determi- 
nation of the nature of these models has therefore required 
much longer simulations to be performed and that no initial 
perturbation was introduced. A summary of these long-term 
simulations is reported in Fig. [161 which shows the evolu- 
tion for models Ul and SI. The upper panel, in particular, 
shows the amplitude of the distortion parameter 77X for model 
SI (continuous line) and Ul (dotted line), respectively, while 
the lower panel shows the evolution of the power in the m=2 
mode; indicated with a dashed line is the average value of the 
TO=4-mode power for model SI, which does not show an ap- 
preciable growth. 

While it is evident that model Ul develops an instability 
over the timescale of the simulation, SI does not, implying 
that the threshold for the onset of the dynamical bar-mode in- 
stability for the sequence under investigation is Pc — 0.255. 
We are aware that it may be argued that model S 1 is also unsta- 
ble but with a much smaller growth rate; assessing in practice 
this would require simulations which are prohibitive with the 
present computational facilities. However, confidence that the 
result obtained is coiTect is also provided by critical-fit analy- 
sis, which will be discussed in the following Section. 



B. Second method: critical fit 

The classical Newtonian study of the bar-mode instabil- 
ity is based on the stability analysis of Maclaurin spheroids 
for an incompressible and self-gravitating fluid in equilib- 
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Model 


p 


52 


ti 


t2 


V 




/b 








(ms) (ms) 


(max) 


(ms) 


(Hz) 


S6 


0.240 


0.04 


3 


9 


0.02 


- 


740t^f 


S5 


0.245 


0.04 


3 


9 


0.02 


- 


705+^f 


S4 


0.250 


0.04 


3 


9 


0.03 


- 




S3 


0.252 


0.04 


3 


9 


0.04 


- 


611+*^ 


S2 


0.253 


0.04 


3 


9 


0.05 


- 




SI 


0.254 


0.04 


3 


9 


* 


9.71 




Ul 


0.255 


0.04 


3 


9 


* 


5.26 


56711 


Si 


0.254 


0.0 


45 


63 


0.02 




= 00+209 

oay_209 


Ul 


0.255 


0.0 


45 


63 


* 


22.1 


588±^« 



TABLE III: The Table reports for some models near the instability 
threshold: the ratio l3 of the rotational kinetic energy to the gravi- 
tational binding energy; the initial m=2-mode perturbation S2; the 
bounds ti and t2 of the considered time interval; the maximum 
distortion 77; the growth rate and the frequency of the bar- 
mode during the initial part of the instability. For these models, 
Ax/AIq — 0.5 and 7r-symmetry is not used. The growth rate 
and the frequency of the bar-mode are obtained making a least- 
square fit of rj+{t) (eg 14.5b between ti and t2. The * indicates that r] 
did not reach a maximum before the simulation was stopped. 



rium [1]. We recall that by considering the linearized form 
of the second-order virial equation, which governs the small 
oscillations around the equilibrium configuration, it is possi- 
ble to show that the toroidal m=2 perturbations have complex 
eigenvalues given by (we here use Chandrasekhar's notation 
of ref. IdJ) 



(7.1) 



where e is the eccentricity of the Maclaurin spheroids relative 
to an incompressible Newtonian star (e = for a spherical 
star) and 



Bn = 



3 e - 5 + 2 e''' + %/l - (4 - 3) arcsin(e) 



4e5 



2 6 (e2 - 1) 2 (3 - 2 e^) Vl - arcsin(e) 

(7.2) 

As customary in perturbative studies, the real part of the 
eigenfrequency a represents the characteristic frequency of 
the small perturbation, while its imaginary part measures the 
exponential growth time of the m=2 mode and is nonzero if 
4Bii(e) - r22(e) < 0, with the threshold for the instabil- 
ity being marked by the value of the eccentricity for which 
4Bii(e) — f2^(e) ~ 0. In the case of an incompressible New- 
tonian star, the eccentricity e and the instability parameter /3 
are related as 



3\/l~^ 



2e^ 2e arcsin(e) 



(7.3) 



Equation ( 17. It is most conveniently rewritten in terms of the 
instability parameter (3 as 
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ti 


t2 


V 




/b 






(ms) (ms) 


(max) 


(ms) 


(Hz) 


U2 


0.2581 


16.9 


22.4 


0.3734 




552+3Q 


U3 


0.2595 


19.9 


24.2 


0.4241 


2.678tlil 




U4 


0.2621 


15.3 


18.3 


0.5496 


1.854tlf, 




U5 


0.2631 


16.2 


19.0 


0.5788 


1.748l»;12o 


538t«3 


U6 


0.2651 


14.5 


17.1 


0.6305 


1.574tni 


528tl 


U7 


0.2671 


14.2 


16.4 


0.6694 


1.408tH« 




U8 


0.2686 


12.2 


14.3 


0.7027 




518tl 


U9 


0.2701 


13.2 


15.2 


0.7223 


1.269lH^ 


512t5 


UIO 


0.2721 


13.7 


15.6 


0.7482 


1 1 84+0-05 


503lf 


Ull 


0.2743 


12.9 


14.7 


0.7749 


r.iro_o.Ol 


493lti 


U12 


0.2761 


12.0 


13.7 


0.7999 


1.0661^^3 


486t| 


U13 


0.2812 


11.2 


12.7 


0.8551 


0.952lH^ 


453t| 



TABLE IV: Same quantities as in Tablelllll but referring to the mod- 
els evolved with a grid size of Ax/Mq = 0.625, with 7r-symmetry 
and no perturbation. The growth rate Tg and the frequency of the 
bar-mode are obtained making a least-square fit of 77+ (t) (eq |4.5| l be- 
tween ti and t2. The interval [ti,t2] is here determined, differently 
from the previous tables, as the one in which rj{t) is between 5% and 
the 25% of its first maximum. These models were used to determine 
the threshold. 



where 



1/7 



-4Bn(/3) + r!'(/3) ■ 



(7.5) 



With this in mind it is natural to ask how accurate is 
the Newtonian description of an incompressible Maclaurin 
spheroid in describing the nonlinear dynamics of a relativistic 
differentially rotating star Of course also in full General Rel- 
ativity the transition between stable and dynamical unstable 
bar-modes will be described by a change from real to imagi- 
nary of the m=2-mode eigenfrequency and it is therefore rea- 
sonable to expect that the frequencies and the growth times 
depend on f3 as 



2tt 



fc + /i') W - Pc) + fP W - Pc)' , (7.6) 



(7.7) 



Using this ansatz, we have performed a series of simulations 
for models U2-U13 using a 7r-symmetry and a grid resolu- 
tion of /S.x/Mq = 0.625, through which we have computed 
the values for and by means of a nonlinear least-square 
fit to the trial form of eq. (I4.5l l. Making use of these results, 
which are collected in Table |IV] we have then again used a 
least-square fit of the values and to the expected (3 de- 
pendence of eqs. ( 17. 6t and ( 17. 7t to obtain the unknown coeffi- 
cients /c, /i^'', k and Pc- 



/c = 554 Hz, /< 
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-1668 Hz, f)^ 



(2)_ 



-85635 Hz, 
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/3c = 0.2554, fc = 0.153 ms 
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FIG. 17: Left panel: critical diagram as constructed with the frequencies and growth times relative to the unperturbed models of Table [W] 
(triangles) and to the perturbed models of Table IIIII (squares). The continuous lines represent the two fitted curves for and r(/9), while 
the dotted line the corresponding extrapolations below the threshold. Right panel: same data as in the left panel but magnified around the 
critical threshold and expressed in terms of 1/r^ to highlight the very good fit. 



In the left panel of Fig. [17] we summarize the results of this 
fit for the models U2-U 1 3 and show the corresponding error 
bars. In particular, we denote with triangles the unperturbed 
models reported in Table II VI while with squares the perturbed 
models in Tablelllll In addition, the continuous lines represent 
the two fitted curves for ri(/3) and t(/3), while the dotted line 
refers to the corresponding extrapolations below the threshold. 
We note that the error bars in Fig.[T7]are computed in different 
ways for the growth rates and the frequencies. In the first case 
they are computed as the difference between the minimum 
and maximum values of d\og{rj{t)) / dt in the time intervals 
between ti and ^2 reported in Table |IV] In the second case, 
instead, the error bars are determined using the minimum and 
maximum values over the time intervals between ti and t2 
reported in Tables [III] and |IV] of the pattern speeds extracted 
from the collective phase 02 (i) of eq. (I4.8l l. 

A number of comments are worth making. Firstly, it is clear 
that the Newtonian description of the instability in terms of 
incompressible Maclaurin spheroids is surprisingly accurate 
also in full General Relativity and for differentially rotating 
stars. It is especially so for models which are far from the 
instability threshold, for which the error bars are very small 
and below 5%. Secondly, as the models approach the critical 
threshold from above (i.e, unstable models), the growth times 
become increasingly large and the numerical errors increas- 
ingly more important for the resolutions we could use. Yet, 
even in this regime the perturbative predictions are accurate to 
better than 15%. Finally, as the models approach the critical 
threshold from below {i.e, stable models), the frequencies of 
the oscillations triggered by the use of an initial perturbation 
are similar to the spin frequencies, making the determination 
of the eigenfrequencies increasingly difficult and hence inac- 



curate (c/ the large error bars for models denoted by squares 
in the left panel of Fig. [TtI i. 

Finally we note that the independent determination of the 
threshold for the instability provided by the fit to eq. ( |7.7| i (i.e. 
(3c = 0.2554) is in surprisingly good agreement with the one 
obtained in the previous Section, i.e. Pc — 0.255, confirm- 
ing the accuracy of the dynamical determination and suggest- 
ing that model SI is indeed stable. The very good fit of the 
data is shown in the right panel of Fig. [iTl which reports the 
same data as in the left panel but magnified around the critical 
threshold and shown in terms of 1 /r^ to highlight the func- 
tional dependence as expressed by eq. ( I7.7l i. 



VIII. GRAVITATIONAL-WAVE EMISSION AND 
DETECTION OF THE UNSTABLE MODELS 

One of the goals of this study is to assess whether a dy- 
namical bar-mode instability triggered in a neutron-star model 
could be a good source of gravitational radiation for the detec- 
tors presently collecting data (LIGO, GEO), in the final stages 
of construction (Virgo) or in the planning phase (Dual) ll42ll . 
The use of uniform grids in these calculations has prevented 
us from placing the outer boundaries at distances sufficiently 
large to allow a gauge-invariant extraction of the gravitational 
waves ll43l l44l |45|] . However, a reasonable measure of the 
amplitude of the expected gravitational radiation and of the 
consequent SNR is still possible by making use of the stan- 
dard quadrupole formula |46]. A number of comparative tests 
have been carried out among the various ways in which the 
gravitational-wave amplitudes can be estimated (see, for in- 
stance, refs. ll23ll47lR8il ) and we expect the error in this case 
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to be 0{M/R) and thus of a few percent only. 

In this approximation, the observed waveform and ampH- 
tude for the two polarizations measured by an observer situ- 
ated at infinity along the 2;-axis are 



r^'it') - iyy{t') 



ry{t') 



(8.1) 



where t' — t — r is the retarded time. The total gravitational 
wave luminosity is then given by 



E 



1 



while the angular momentum loss is 
2 



J,. 



:{^ijki''^ ilk) 



(8.2) 



(8.3) 



We note that the definition of the quadrupole tensor is not 
unique as p and D coincide in a Newtonian approximation and 
either of them could be used. Following ref. il9ll . we adopt the 
definition expressed in eq. ( 14.1b . as this allows us to exploit 
the conservation equation for D to perform analytically the 
first time derivative of the quadrupole tensor and obtain that 

= I d\D [x\av^ + 13^) + x\av' + [3')] , (8.4) 



from which we compute the second and third / time- 
derivatives using first-order finite differencing. This approach 
not only is simpler, but it is indeed the only possible one, since 
our second-order-accurate scheme would not allow for an ac- 
curate direct calculation of / (see discussion in Appendix 
A of ref. iH). 

In Fig. [18] the we show the waveforms for the "cross" po- 
larization computed in this way, as measured along the z-axis 
for models U3, Ull, U13 and D2 (We recall that the "cross" 
and "plus" waveforms for a stationary and rotating bar differ 
only by a phase, since the emission is circularly polarized). 
Clearly, the gravitational-wave signal mimics (modulo a fac- 
tor of 2 in the frequency) the dynamical behaviour of the bar 
deformation, with signals that are longer for models near the 
stability threshold {i.e. model U3), and with amphtudes which 
increase for the more compact model D2. 

A convenient measure of the strength of the signal can be 
given in terms of the root-sum-square amplitude of, say, the 
plus polarization 



4-00 



dt h%{t) 



1/2 



(8.5) 



since /ijss has the same units as the strain-noise amplitude of 
the detectors and it is therefore possible to obtain a rough es- 
timate of the SNR simply dividing it by the strain-noise am- 
plitude at the frequency where the signal is the strongest. In 
Table [V] we report the values of h^^s for a signal coming from 
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FIG. 18: Gravitational-wave signals along the 2:-axis for the "cross" 
polarization as computed in the Newtonian quadrupole approxima- 
tion for models U3, Ull, U13 and D2. 



Model 


AM/M AJ/J 

(10^^) (10^^) 


^rss 
(10"^") 


SNR SNR SNR SNR 
V L AL Dual 


U3 


7.75 


3.82 


1.39 


117 


82.5 1590 19 


Ull 


9.48 


4.84 


1.74 


155 


119 2040 23 


U13 


6.23 


3.18 


1.51 


138 


114 1780 30 


D2 


52.6 


10.2 


2.50 


136 


79.8 2270 892 



TABLE V: List of the representative gravitational-wave quantities 
computed in the Newtonian quadrupole approximation for models 
U3, Ull, U13 and D2. From the left the different columns report: 
the fractional amounts of the energy and angular momentum carried 
by the gravitational radiation {AM/M and A J/ J, respectively), the 
root-sum-square of h+ for a source at 10 kpc, and the SNRs for 
Virgo, LIGO, Advanced LIGO and Dual. 



a source at 10 kpc, together with the SNR computed assuming 
the optimal use of match-filtering techniques and given by 




+ 00 



dlog/ 



\h+{f)?f 

SnU) 



(8.6) 



where /i+ (/) is the Fourier transform of (t) and S'„ (/) is 
the designed sensitivity of either Virgo (V), LIGO (L), Ad- 
vanced LIGO (AL) or of the planned resonant detector Dual. 
Of course, the SNR is inversely proportional to the distance 
from the source and the values reported in Table [V] indicate 
that, while an instability developing in a rapidly rotating star 
in our Galaxy would yield an extremely strong signal, the 
SNR is still expected to be 0(1) also for a source as far as 
about 10 Mpc if measured by Advanced LIGO. A different 
representation of the results summarised in TablelVlis offered 
in Fig. [19] which shows a spectral comparison between the 
designed sensitivity hrms (/) = yS^Jj) of Virgo, LIGO and 
Advanced LIGO and the power spectrum \h+{J)\f^ of the 
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FIG. 19: Comparison between \h{f)\f^'''^ for models U3 and U13 at lOkpc and the square-root of the power spectrum of the noise of Virgo 
(dashed hne), LIGO (dotted line), Advanced LIGO (dotted-line) and the planned resonant detector Dual (dot-dashed line). Note the significant 
difference in the power spectrum of the two signals, with the one relative to model U3 having a larger and narrower peak at about 600 Hz, 
produced by the more persistent bar deformation. 



expected signals for models U3 and U13. 

It is interesting to note the significant differences in the 
power spectrum of the two signals, with the one relative to 
model U3 having a larger and narrower peak as a result of a 
more persistent bar. Yet, the SNR for model U3 is signifi- 
cantly smaller than the one for both models Ul 1 and U13 and 
a number of different factors are behind this somewhat sur- 
prising result. Firstly, while the bar is indeed more persistent 
for U3, the amplitude of the bar distortion is larger for models 
Ull and U13, thus larger is the corresponding gravitational- 
wave signal. Secondly, the frequency of the power-spectrum 
maximum is smaller for models Ull and U13 and thus bet- 
ter fitting the sensitivities of the detectors. Finally, the ide- 
alised assumption that match-filtering techniques can be used 
at all frequencies implies that all of the power spectrum con- 
tributes to the final SNR and the large wings of the spectra 
of models Ull and U13 can significantly increase the SNR. 
Indeed, when evaluating eq. (18.6b in a narrow window of 100 
Hz around the peak, the SNR of the different models is com- 
parable. 

Of course, a strong SNR is just a necessary condition for the 
detection and the possibility of measuring the gravitational- 
wave signal from this process will depend significantly on 
its event rate, which is still largely uncertain. In the case in 
which the instability develops in a hot protoneutron star re- 
sulting from the collapse of a stellar core, the event rate is 
strictly related to the supernova event rate, which is of 1 or 2 
supernova per century per galaxy. About 60% of the remnants 
of the explosion should be neutron stars, but the requirement 
of rapid rotation in the progenitor makes the event rate of dy- 
namical instabilities considerably lower |50], although such 
prospect may be more optimistic after recent studies L5lll . On 



the other hand, in case the bar is produced as a result of a 
binary merger of neutron stars, the most optimistic scenarios 
suggest that such mergers may occur approximately once per 
year within a distance of about 50 Mpc |52]. Finally, the event 
rate of the classical scenario in which the instability is trig- 
gered in an old neutron star spun up by accretion in a binary 
system, still remains difficult to quantify. 

IX. CONCLUSIONS 

We have presented accurate simulations of the dynamical 
bar-mode instability in full General Relativity. An impor- 
tant motivation behind this work is the need to go beyond the 
standard phenomenological discussion of the instability and 
to find answers to important open questions about its nonlin- 
ear dynamics. Among such open problems there are, for in- 
stance, the determination of the role of the initial perturbation 
or of the symmetry conditions, or the influence on the dynam- 
ics of the value of the parameter /3, or, most importantly, the 
determination of the timescale of the persistence of the bar 
deformation once this is fully developed. Clearly, this latter 
question is a very pressing one in gravitational-wave astron- 
omy, as it bears important consequences on the detectability 
of the whole process. 

In order to provide answers to these questions we have 
explored the onset and development of the instability for a 
large number of initial stellar models. These have been calcu- 
lated as stationary equilibrium solutions for axisymmetric and 
rapidly rotating relativistic stars in polar coordinates. More 
specifically, the evolved models represent relativistic poly- 
tropes with adiabatic index F = 2 and are members of a se- 
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quence having a constant amount of differential rotation with 
A ~ 1 and a constant rest mass of Mq ~ 1.51 Mq. 

The simulations have been carried out with the general- 
relativistic hydrodynamics code Whisky, in which the hy- 
drodynamics equations are written as finite differences on 
a Cartesian grid and solved using HRSC schemes ||25l l26ll . 
The Einstein, equations, on the other hand, have been 
solved within the conformal traceless formulation imple- 
mented within the Cactus computational toolkit |30']- 

The main results of our analysis can be summarised as fol- 
lows: i) An initial m=l or m=2-mode perturbation can play 
a role in determining the duration of the bar-mode deforma- 
tion, but not the growth time of the instability; the only ex- 
ception to this is represented by models near the threshold; ii) 
For moderately overcritical models the use of a 7r-symmetry 
can radically change the dynamics and extend considerably 
the persistence of the bar; this ceases to be true for largely 
overcritical models; Hi) The persistence of the bar is strongly 
dependent on the degree of overcriticality and is generically 
of the order of the dynamical timescale; iv) Generic nonlinear 
mode-coupling effects (especially between the m=l and the 
m=2 mode) appear during the development of the instability 
and these can severely limit the persistence of the bar defor- 
mation and eventually suppress the bar deformation; v) The 
dynamics of largely overcritical models (i.e. with /3 ^ /3c) are 
fully determined by the excess of rotational energy and the bar 
deformation is very rapidly suppressed through the conversion 
of kinetic energy into internal one. Interestingly, a similar dy- 
namics for the odd and even modes has been observed also in 
the case of the low-/3 instability in recent Newtonian simula- 
tions |53]. In this case, however, the growth rate of the m=l 



mode is much smaller and hence the persistence of the bar 
longer 

Finally, we have considered whether the classical Newto- 
nian stability analysis of Maclaurin spheroids for incompress- 
ible self-gravitating fluids is accurate also for differentially ro- 
tating and relativistic stars. Overall, we have found the pertur- 
bative predictions to be surprisingly accurate in determining 
the threshold for the instability as well as the complex eigen- 
frequencies for the unstable models. 

While the features of the bar-mode instability discussed 
here are expected to be rather generic, they have been deduced 
from the analysis of a small region of the space of parame- 
ters. Work is now in progress to extend the present analysis 
by considering stellar models with different and larger com- 
pactnesses and that are regulated by more realistic EOSs. 
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